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struction procedure based again on numerical fluxes rather than cell- averages. We also consider two 
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plementation of the ENO schemes. Two methods of sharpening contact discontinuities-the subcell 
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I. INTRODUCTION 


In this paper we extend the construction we began in [12] of efficient imple- 
mentation of ENO (essentially non-oscillatory) schemes approximating systems of 
hyperbolic conservation laws of the type: 
d 

u t + 53 = 0 (or = g(u,x, i), a forcing term) (1.1a) 

i=i 

u(x,0)=u°(x) (1.1b) 

where u = (uj, . . . ,u m ) T , x = (s 1 ,® 2 ,...,® 4 ), and, for real £ = (£i,...,£<*), 
the combination £»§{[■ always has m real eigenvalues and a complete set of 

eigenvectors. We use the notations Xj = J-Ax, t n = nA£, and use Uj to denote 
the computed approximation to the exact solution u(xj,t) of (1.1). We use bold 
face letters for vectors and plain letters for scalars. 

ENO schemes, constructed by Harten, Osher, Engquist, and Chakravarthy [2], 
[3], [4], [5], use a local adaptive stencil to obtain information automatically from 
regions of smoothness when the solution develops discontinuities. As a result, ap- 
proximations using these methods can obtain uniformly high order accuracy right 
up to discontinuities, while keeping a sharp, essentially non-oscillatory shock tran- 
sition. The original ENO schemes in [2], [3], [4], [5] used a cell-average framework 
which involved a reconstruction procedure to recover accurate point values from 
cell averages, and a Lax-Wendroff procedure (replacing time derivatives by space 
derivatives, using the PDE) for the time discretization. This can become a bit 
complicated for multi-dimensional problems [2]. Harten is currently investigating 
efficient and more local methods [7]. For ease of implementation we constructed [12] 
ENO schemes applying the adaptive stencil idea to the numerical fluxes and using 
a TVD Runge-Kutta type time discretizations. The ENO schemes in [12] skip the 
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reconstruction step and the Lax-WendrofF time discretization procedure, hence they 
are simpler to program, especially for multi-space dimensional problems. We hope 
that [12] and this paper will encourage potential users to apply these high order 
methods to real physical computations for systems of conservation laws in several 
space dimensions. 

We shall use the same notation as in [12]: A± are the usual difference operators 
A ±ay = ±(oy±i — ay); the equation (l.la) is sometimes written in an abstract form: 

u t = jC(u) (1.2) 

The ENO spatial operator L(u) is supposed to approximate £(u) to r-th order: 

L(u) = £(u) + 0(h r ) (1.3) 

for smooth u, where h is the maximum mesh size, and the Euler forward version 


w = T(u) = (I + AtL) (u) 


(1.4) 


is assumed to be total variation stable, for scalar, ID nonlinear problems, under a 
suitable CFL restriction 


A = < A 0 

Ax 


(1.5) 


where Ao is usually inversely proportional to max |/ , (u)|. At present this stability 
cannot be proven for unmodified third order or higher ENO schemes. There is, 
however, strong theoretical and numerical evidence to indicate that the methods 
are indeed TV stable [3], [4], [5], [12]. 


An r-th order TVD Runge-Kutta time discretization is then applied: 

u (,) = J2 [ a »'* u(fe) + ft*AtZr(iiW], * = 1, . . . , f 

k = 0 


(1.6a) 
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u (°) _ u n j u (?) _ u "+l (1.6b) 

with 

TV(u n+1 ) < 2V(r(u n )) (1.7a) 

under the CFL restriction: 

At 

A = — < c r • A 0 . (1.7b) 

We call c r the CFL coefficient. We constructed schemes [12] for r < 4 with f — r 
and C2 = C3 = 1, c\ = |. For r = 5 we needed f = 6 and C5 = For r > 4 we 
also needed L which approximate £, in an “adjoint” equation u t = —C( u) because 
some of the /?,-* become negative. The details are presented in [12]. Some of the 
schemes used in Section 5 are listed in Table 1 at the end of this paper. 

We shall always use conservative schemes, i.e. for the scalar ID problem (1.1) 
we write: 

= (i-N 

for a consistent numerical flux 

= /(«;-*•••.«;+*); /(u,...,u) = /(u). (i-8b) 

We proved in [12] the existence of constants ai, 04 , • • • such that 

[ S ! a** 

fj+k = f i+7 + £ a 2fcAx 2i: (^/) y+ i + 0(Az r+1 ) (1.9) 

k—1 

guarantees r-th order accuracy (1.3). For example, 02 = — a 4 — 5^, • • • 

The use of the a 2 k may seem a bit unnatural. We need to evaluate all the even 
derivatives of the interpolating polynomials (except the last one if r is even), which 
increases the computational cost for large r. In Section 2 of this paper we provide 
a simplified (but equivalent) version of (1.9). We also consider in Section 2 two 
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improvements labelled ENO-LLF (local Lax-Friedrichs) and ENO-Roe, which yield 
sharper shock transitions, better overall accuracy, and lower computational costs 
than the methods described in [12]. 

Shocks have a self-sharpening mechanism due to converging characteristics. 
This is not true for contact discontinuities, which are usually smeared in shock 
capturing calculations. Recently, Harten [6] introduced the notion of subcell reso- 
lution in order to sharpen contact discontinuities. He did this in the context of cell 
average based ENO schemes. Although in [6] the fact that cell averages contain 
information leading to the location of the shock is strongly used, this information 
is actually contained in any conservative scheme (1.8). In Section 3 we translate 
Harten’s subcell resolution idea to the ENO schemes using fluxes and Runge-Kutta 
techniques. The result is simpler than [6] but less accurate in theoretical resolution. 
However the numerical results (in Section 5) are quite close to those presented in 
[6]. The main drawback of our extension of Harten’s approach is that it is not clear 
how to generalize it effectively to several space dimensions. A naive generalization 
gave unsatisfactory results in our 2D numerical tests. Currently Harten is consid- 
ering truly 2D subcell resolution methods [7]. Also recently Yang [15] introduced a 
simple artificial compression technique applied to cell average based ENO schemes. 
He demonstrated that the technique preserves the TVD (or ENO) properties when 
applied to a TVD (or ENO) schemes. We translate his result to our present frame- 
work in Section 3. The numerical results presented in Section 5 indicate that this 
adaptation works well both for ID and 2D problems. 

In Section 4 we extend the present implementation to nonlinear systems and to 
multi-space dimensions. Section 5 contains numerical experiments, including scalar, 
ID nonconvex Riemann problems; examples of sharpening contact discontinuities; 



2D Riemann problems for Burgers’ equation; ID and 2D Euler equations of gas dy- 
namics including a compressible shock-turbulence interaction problem. Concluding 
remarks are presented in Section 6. 


n. SIMPLIFICATIONS AND IMPROVEMENTS 

We begin with the ID scalar nonlinear problem, i.e. d = m = 1 in (1.1). 

To simplify (1-9) we use the following elementary result: 

LEMMA 2.1. If a function h(x) satisfies 

r x i 

/(«(*)) = (~ 1 ) 

then /(«(*)), = M " 1 •¥>-><--¥> □ 

The proof is trivial. 

Lemma 2.1 tells us that the numerical flux should approximate k(xy + i) 
to a high order. It is not easy to obtain h(x) directly from (2.1). However, the 
“reconstruction via primitive function” technique in [4] can be applied to obtain 
the primitive function of h(x) 

H(x) = f h(t)dt (2.2) 

J — OO 

at Xy+ 1 . by 

H(x j+i ) = [’** h(()i( = £ h(()d( = it £ /(«») (2.3) 

00 k=—oa X k— y fc= — OO 


Notice that the lower limit — oo is irrelevant. It can be changed to any fixed 
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Given the values of H(xj + L), we can construct interpolating polynomials in 
an ENO fashion, i.e. by obtaining a locally “smoothest” stencil starting from one 
or two points, then adding one point to the stencil at each stage by comparing two 
divided differences and choosing the one which is smaller in absolute value, (see 
[4]), and then use its derivative at Xy + 1 . as the numerical flux fjH- We remark 
here that by (2.3), we have + ^ = /(tty), hence the ( k + l) 3t order 

divided differences of H can be easily obtained from the fc-th divided differences of 
/. Since we shall never need the 0-th order divided differences, we do not perform 
the su mm ation in (2.3). We simply use the divided difference tables for /. 

Lemma 2.1 also explains why the ENO schemes constructed here and in [12] 
are equivalent (for linear equations in one space dimension only) to the cell-average 
ENO schemes constructed in [4] via the primitive function for u. 

The starting point in the choice of stencil process is very important. “Up- 
winding” is achieved by this initial choice, and this also is crucial for the evident 
stability of these methods. We have experimented with using x/+i as the starting 

A 

point in the choice of stencil process for evaluating / J+ i., to solve Burgers’ equation 
(i.e. upwinding in the wrong direction initially) with positive f'[u). We observed 
instability even if the stencil is chosen in an ENO fashion up to third order. In [12] 
/ is decomposed into / = /+ + /“ with (/ + )' > 0, (/“)' < 0, then different start- 
ing points are assigned to /* according to the direction of the wind. This simple 
procedure may smear shocks and affect the overall accuracy (although the smearing 
is very mild comparing with lower order non-ENO schemes). It was also pointed 
out in [12], that instead of decompositing / into / + and f~ we could also use any 
E-flux hj + 1 . [10], with dff + ^ = /y +1 —h }+ ^, df~ + ^ = \ ~ fj to replace the first 

(undivided) differences f± +1 — we shall later refer to this procedure as ENO 
schemes with building blocks. However that procedure requires smoothness 



-7- 


of h :H at sonic points to keep the high accuracy there. It is known that most in- 
fluxes are not smooth at sonic points, they are typically only Lipschitz continuous 
or perhaps C 1 there). Although wc may modify any £-fiux to make it smooth at 
sonic points, this procedure is too complicated to be computationally pleasing. 


We now use a different approach to overcome this difficulty. We choose the 
first point in the stencil according to the local sign of /'(u) at x J+ i,. We may use, 
e.g., the “Roe” speed 


a 3+k = 


/( u j+i) ~ /(« j) 

Uj+i - Uj 


(2.4) 


to determine the sign of /'(u J+ i,). We then have the following algorithm: 


Algorithm 2.1. (ENO-Roe). 


(1) Compute the divided difference table of /*, and identify 


H[x t _^x t ^} = f[u{x t ) } 


(2.5a) 


H[x e _±,x l+ ±,...,x t+k+ ±\ 


k 

k + 1 


f[u[x t ) 


,...,u(x i+k )], k = 1,2, . . . , r 


(2.5b) 


(2) If a,- , jl = > 0, then 

v > J-t-3 «i+I-»y — ’ 


mm 


(2.6a) 


else 


=7 + 1 
K mm J ^ 


(2.6b) 


(3) Q (1) (x) = IT[s fc <i> + ^](x 


X,(i) i) 
''min - ? 
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(4) Inductively, if and Q^ -1 )(x) are both defined, then let 

“ w = 

i ( ‘> = Elz k «-„_ 1 _ i 


and 


(i) if ja|(^ > |6^|, then 


c<‘> = »«, 41 = 4l‘> - 1 


(2.7a) 

(2.7b) 


(2.8a) 


otherwise 


C W = a W, fc (£) = 

’ "min min 


(2.8b) 


(ii) form 




gW(s) = ^" 1) ( l ) + « w E[ (*“**-*) 




(5) $,+*(*) 


(2.9) 


(6) we then take 


fl+i = ii Q i+>iW I—, 


j+i’ 


( 2 . 10 ) 

□ 


REMARK 2.1. (a) Notice that unlike the ENO scheme which uses the Roe flux 

A 

as a building block, Algorithm 2.1, gives a uniformly high order flux f ]+ i, , even at 
sonic points. The scheme described above will be (r + l) 3t order accurate except 
perhaps at isolated zeros of derivatives of the flux f(u(x)) where it may degenerate 
to r-th order. 
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(b) Since x = x J+ j. is a node of Q } +t in (2.10), the evaluation of the derivative 
in (2.10) costs very little. □ 

Algorithm 2.1 yields a scheme admitting a stationary entropy- violating expan- 
sion shock since it is based on the first order Roe scheme. However, the “entropy 
fix” in this framework turns out to be very simple. Before stating this remedy, let 
us present another algorithm: 

Algorithm 2.2. (ENO-LLF) 

(1) Compute the divided difference table of / and u, and identify 

s± [ x t-^x t+ ^\ = |(/[u(x*)] ± a y+ £u[x*]), (2.11a) 

l = / — r, . . . , j + r for H + 
and t = j — r + 1, . . . , j + r + 1 for H~\ 

H ± [x l _± i x t +±,...,x t+k+ ).\ (2.11b) 

i = j — r, . . . , j -+• r — k for H + 

and L = j — r + 1,.. . ,/ + r — k + 1 for H~ , k = l,2, ...,r 

where 

otj+L = max |/'(u)| (2.12) 

(2) For H + , k^ n = j, then repeat steps (3)-(4) in Algorithm 2.1 to get 

Qi r+1) (x). Let Qf +i (x) = <?<T +1 >(x): 
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(3) For H~, = j + 1, then repeat steps (3)-(4) in Algorithm 2.1 to get 

OL r+1) (x). Let Qj +i = <3i r+1) (x); 

(4) We then take 

+ < 2 ' 13 > 

□ 

REMARK 2.2. (a) When / is convex, f" > 0, (2.12) simplifies to 

«y+i- =max(|/'(uy)|, |/ # (u y +i)|) (2.14) 

(b) It can be verified that the local Lax-Friedrichs flux 

h)+\ = \[f{uj) + /(« y+i) ~ a y+ j( u J+i “ u j)1 ( 2 - 15 ) 

where a y+ 1 . is defined by (2.12), is a monotone flux (see [1] for definitions). We omit 

the details. However Algorithm 2.2 is not equivalent to the ENO scheme obtained 

by using as a building block [12]. The advantage of Algorithm 2.2 is again its 
J"r y 

uniformly high order of accuracy; 

(c) In terms of cost, Algorithm 2.2 is equivalent to ENO schemes using in 
[12] (we shall denote it by ENO-LF), i.e. evaluating two divided difference tables, 
while Algorithm 2.1 is half as expensive because it only needs one divided difference 
table; 

(d) Roughly speaking the viscosity of ENO-LF is larger than that of ENO-LLF, 
which is in turn larger than that of ENO-Roe. Hence in going from ENO-LF to 
ENO-LLF to ENO-Roe we can expect increasingly less shock smearing and better 
overall accuracy; this is verified by the numerical experiments of Section 5. 
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(e) As was true for the schemes described in [12], when r = 1 Algorithm 2.1 
yields a familiar TVD second order accurate scheme [9], when the minimum of 
absolute value function is changed to the minmod function: 

{ s • min |a,-|, if s = sign (ai) = • • • = sign (a n ) 

1 -‘- n (2.16) 
0, otherwise 

We omit the details of the derivation here. □ 


Algorithm2.Z (ENO-RF). Let fj+h be defined by (2.10) if /'(u) does not change 
sign between uy and uy+i; otherwise, let be defined by (2.13). □ 

REMARK 2.3. If / is convex, /"(u) > 0, we can use (2.10) unless /'(uy) < 
0 < /'(uy + i), i.e. we use (2.13) only in the cells which contain “expansion shocks”. 
This is the case for, e.g. the Euler equations of fluid gas dynamics. See Section 5. 
□ 

Notice that the cost of Algorithm 2.3 is very close to that of Algorithm 2.1 
because sonic points are isolated and the divided differences of u need only be 
computed locally near sonic points. 


We may associate Algorithm 2.3 with the following E-flux. 



f /( u j). 

/(“/+ 1 )> 

*isj. 


if /'(u) > 0 between uy and uy + i 
if f'(u) < 0 between uy and uy+i 
otherwise 


(2.17) 


however, as mentioned above this algorithm is not equivalent to the ENO scheme 
obtained by using as a building block. 


Numerically we observe that Algorithm 2.3 always yields the correct entropy 


solution even for nonconvex /. See Section 5 for some examples. 



E3. SHARPENING OF CONTACT DISCONTINUITIES 

We first consider the ID, scalar, linear version of (1.1), i.e. d = m = 1, 

/(u) = au. Discontinuous solutions to this problem (called contact discontinuities) 
are usually smeared more severely than shocks, when computed by shock capturing 
schemes. 

Recently Harten [6] introduced the concept of subcell resolution which he used 
to sharpen contact discontinuities. This notion uses the observation that the cell 
averages carry information about shock locations. In fact, this information is also 
contained in any conservative scheme (1.8), whether cell average or numerical flux 
based. We can thus translate Harten’s subcell resolution techniques to our ENO 
schemes which use fluxes and Runge-Kutta time discretizations. Notice that for a 
linear, constant coefficient problem /(u) = au, ENO-LF, ENO-LLF and ENO-Roe 
are identical. We write out the following algorithm for the case a > 0. 

Algorithm 3.1. For a > 0, at the beginning of every Runge-Kutta cycle: 

(1) Let Sj = |m(A + u" , A_u" )j where m is the minmod function (2.16). Define 
the “critical intervals” (intervals containing discontinuities) Jy = (xy_ i., x ;+ i.) 

by Sj > Sy+i, Sj > Sj—i] 

(2) For any “critical interval” Jy, let Oj = - or more accurately let Oj be 

the solution in [0,1] of the quadratic equation 

(uy+ 1 +tty— i-Uy + 2 -Uy— 2 )0y+(3uy +2 +3uy_ 1 -5uy +1 -Uy_2)0y+(4uy +1 -2uy + 2— 2uy) 

we use Xy_i. + 0yAx as an approximation to the discontinuity location inside 
the cell Jy; 



We then perform each step of the Runge-Kutta cycle as follows: 


(3) Let fj+L be defined as usual (Algorithm 2.1) unless Ij or (for the second, 
third, etc. steps in the Runge-Kutta cycle) Ij-\ is a “critical interval”. If Ij is 
a “critical interval” , we define 


fj+k = ( 1 ~ Z])fj+\ + tjfj+l with b = min ( 


(fl) 


(1 — 0j)Ax 
aAt 


,1) (3.1) 


where fjQ and fjQ are again computed by Algorithm 2.1 with the modi- 
fication on step (2) only i.e., the first step in the choice of stencil. We use 
^min = 3 ~ 1 for fj+\ and k£) n = j + 1 for To be safe, for the second, 

third, etc. steps in the Runge-Kutta cycle, we may choose k^ n = j + 2 for 
f[^\ if £j < 1, and, when Ij- i is a “critical interval” and £j_i < 1, we should 

Jt* y 

choose = j + 1 to evaluate fj+±- □ 

REMARK 3.1. (a) The case a < 0 is easily obtained by symmetry. 


(b) The discontinuity detector in step (l) is somewhat different than that of [6]. 
The purpose of our choice is to avoid identifying an interval containing a smooth 
extremum as “critical interval” . 


(c) As in [6], the modification in (3.1) does not affect accuracy if u happens to 
be smooth in Ij. 

(d) The whole procedure here is simpler than in [6], especially if the first 
expression for Oj is used in step (2). The price we pay is that it is less accurate 
in theoretical resolution, because it is based on piecewise constant (first choice of 
Oj) or piecewise linear (second choice of Oj) interpolants. Numerically we observe 
results very close to those found in [6]. See Section 5 for details. 
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(e) The naive generalization to 2D, namely writing u + A tL x u + A tL y u as 
|(u + 2A tL s u) + \{u + 2A tL y u) and applying the algorithm separately to the two 
parts, does not work satisfactorily, according to our numerical tests (see Section 
5). A truly 2D subcell resolution seems to be needed. Currently Harten [7] is 
investigating this within the cell average framework. 


(f) For details of the analysis of Algorithm 3.1, we refer the readers to [6]. Al- 
though the two algorithms are different, we have essentially borrowed the philosophy 
of [6] and applied it to the present case. □ 


Another interesting method we can use to sharpen contact discontinuities is 
an artificial compression method introduced by Yang [15], again in a cell average 
context, i.e. for (MUSCL type schemes). Translated to our current framework, it 
gives the following algorithm (again for a > 0) 


Algorithm 3.2. Let f* + ^ — f J+ i. + c J+ ^ where 


c /+t - 4'+t’ A'-i fj+k’fj- -fj- 1 ! ( 3 - 2 ) 

Here, as in (3.1), f^\ corresponds to = j + 1 in (2.6), step (2) of Algorithm 


?(«) 


(K) 


2.1. m is again the minmod function (2.16). 


□ 


REMARK 3.2. (a) a in (3.2) is a positive parameter. It can be tuned for each 
individual problem to optimize the results. We usually use 2 < a < 5, but have not 
yet found a general rule for determining a. 


(b) The approximation for a < 0 is constructed symmetrically. 


(c) Other versions of (3.2) are possible. For example, we may use 

V P 


OC = do I 


[uy+i - 2 tty + tty-i| 
|uy + i - Uyj + |uy - Uy_! 


(3.3) 
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in (3.2), where a o > 0 and (3 = integer > 0 are two parameters. The heuristic 
reason to use (3.3) is that a is then very small (~ 0(/i^)) except near a discontinuity 
in either u or its derivative, hence over-compression in smooth regions should be 
avoided. Computationally this improves the performance, especially with some 
tuning of (3. We have not yet been able to find a general rule in choosing the 
optimal /?. □ 

Both algorithms in this section are easily generalizable to variable coefficient 
or nonlinear problems. However special caution is needed when one tries to sharpen 
a (nonlinear) shock, to avoid obtaining a nonphysical, entropy condition violating 
solution. In the computation for the Euler equations of compressible gets dynamics 
(Section 5), we only use the compression in the linearly degenerate fields. 

IV. IMPLEMENTATIONS IN MULTI-DIMENSIONS AND SYSTEMS 

A special advantage of ENO schemes using fluxes and Runge-Kutta methods is 
their relative simplicity in multi-dimensions. The algorithms in Section 2 are applied 
to each of the terms fi(u) x . in (l.la), keeping all other variables fixed. The Runge- 
Kutta method (1.6) is then applied. A typical CFL restriction ^ max) /'(it)! < Ao 
will be replaced by At max^iLi 3 ^ 7 l/i( u )l < Ao. 

For nonlinear systems, we simply apply the algorithms in Section 2 in each 
(local) characteristic field. We take a ID system to exemplify this process. Let A ]+ l 
be some “average” Jacobian at Xy + i. Examples include Aj + i = f£| u _£( u . +u , +1 ) 
or, in the case of Euler equations of gas dynamics A J+ i = f£| u _ u <*°<o where 

= R( u /> u /+i) » s the Roe average of uy and uy+i [10]. We denote the 
eigenvalues and left and right eigenvectors of A,- + i by A^,, , r^,, p = 

1 , . . . , m, normalized so that = 6 pq = 1 if p = q = 0 if p # q. For any 
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vector a, 


fl ( P ) =£ (p) . a 

J + T 


is the component of a in the p-th (local) characteristic field, because 


m 



p=l 


Algorithm 2.1 now becomes: 


(4.1) 


(4.2) 


Algorithm 4.1. 


(1) Same as step (1) in algorithm 2.1, changing vectors to bold face letters; 

(2) For each j, compute . . . ,x^ +fc+ i] for £ = j — 7, . . . , j + 1 and k = 

1, 2 , . . . , r, by (4.1). Here 1 < p < m; 

(3) Apply steps (2)-(6) of algorithm 2.1 to using a,- + 1 = A^i , to get ; 

2 J 1 j J j y 

(4) Use (4.2) to get 4+* □ 

We may similarly generalize Algorithms 2.2 and 2.3. Notice that (2.12) be- 
comes, for the p-th field 



max 

ueL(u,-,u J+ ,) 


|A w (o)| 


(4.3) 


where L(uj-,u j+1 ) is some curve (e.g. a straight line) in phase space connecting 
u j and Uy +1 . For Euler equations of gas dynamics, the fields are either genuinely 
nonlinear or linearly degenerate, hence we may use 


= max(|AW(uy)|, (A^>( U y +1 )|), (4.4) 


which is similar to (2.14). 



NUMERICAL RESULTS 


We use the notation ENO —X—(Y)r where X = Roe, LLF, RF or LF refers to 
Algorithms 2.1, 2.2, 2.3 or the /* version in [11]; Y = S or A refers to Algorithms 
3.1 or 3.2; r is the order of the scheme. For the TVD Runge-Kutta methods used 
please see Table 1. We have run most examples using second, third and fourth order 
schemes, but here we usually include only third order results as representatives. 


EXAMPLE 1. This is the Example 4 in [12, Section IV] revisited. We solve 
the scalar ID Riemann problems of (1.1) with two nonconvex functions /. 


/(«) = j (« 2 - D (« 2 - 4 ) 

4u 2 

/(U) = 4u 2 + (1 - u) 2 


(5.1) 

(5.2) 


ENO-RF-3 is used. The results are displayed in Figures 1, 2, and 3. We observe 
considerable improvements in shock transition and overall resolution compared to 
the results shown in Figures 14, 18, and 22 in [12], ENO-RF-3 and ENO-RF-4 were 
also tested on several other nonconvex problems. We always observed convergence 
to the correct entropy solution and a well resolved solution. When applied to 
Burgers’ equation (Example 1 in [12, Section IV]) we observed improvements in 
accuracy and shock transition of ENO-RF over ENO-LF. We omit the details here. 
□ 

EXAMPLE 2. This example uses the two contact-discontinuity-sharpening 
algorithms in Section 3 applied to scalar ID linear problems. We solve the model 
equation 


ut + u z =0 — 1 < x < 1 
u(x, 0) = u°(x), u°(x) periodic with period 2 


(5.3a) 

(5.3b) 
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Notice that in this linear case all the three algorithms ENO-Roe, ENO-LLF 
and ENO-RF in Section 2 are identical. 


Four initial conditions u°(x) are used. The first three are used by Zalesak [16], 
using 100 equally spaced grid points in [-1, 1) with 


« ri, — r < x < 4 

u°(x) = < 5 5 , 800 time steps, CFL = 0.2 

i 0, otherwise 

u°(x) = e -300z3 , 600 time steps, CFL = 0.1 

„ 0 (I) = {(!-(¥*)’)>. W<* 


(5.4a) 

(5.4b) 


10 ,600 time steps, CFL = 0.1 (5.4c) 

otherwise 


(since u°(x) is periodic we only give the definitions in [-1, l)). The last is used by 
Harten et al [4], [6] 


u°(x + 0.5) = / 


— xsin(§7rx 2 ), 

| sin(27rx)|, 

2x = 1 — sin(37rx)/6, 


1 — x < 3 

|x| < | t = 2 and t = 8. (5.4d) 

i<x<l 


The results are displayed in Figures 4-18. We observe apparent improvements 
of both algorithms in all cases. For comparisons we refer the readers to [6], [16]. □ 


EXAMPLE 3. We solve a Riemann problem for the 2D Burgers’ equation 



of the type: 


u(x,t/,0) 


' Ul, 

X > 0, 

y > o 


x < 0, 

y > o 

«3, 

x < 0, 

y < 0 


x > 0, 

y < 0 


(5.5a) 


(5.5b) 


Depending on the orders of the u,’s, there are eight essentially different solu- 
tion types. See [13] for details. We used ENO-RF and observed convergence to 
the correct entropy solution with good resolution for all cases. Some results are 
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displayed in Figures 19-28. For one case we present the results of ENO-RF-3 for 
20 x 20 and 80 x 80 grid points, and the result using first order Engquist-Osher 
scheme for 80 x 80 grid points. We observe that ENO-RF-3 using 20 x 20 points has 
roughly the same resolution as does the first order montone scheme with 80 x 80 
points. For the remaining seven cases we only show the results of ENO-RF-3 with 
80 x 80 points. 

We also remark here that Yang’s artificial compression Algorithm 3.2 improves 
the resolution in this example according to our numerical experiments. Since we 
are mainly interested in applying Algorithm 3.1 and 3.2 to linear problems, we omit 
the details here. □ 

EXAMPLE 4. This example is suggested by Professor D. Gottlieb (private 
communication). It simulates a 2D boundary layer problem. We solve the boundary 
value problem 

u t + (-r-)x + u y = 0 0<x<2 7T, 0 < y < 1 (5.6a) 

St 

u(x, 0,£) = a + /Jsinx, u is periodic in x with period 2n (5.6a) 

to steady state. In steady state it resembles the ID Burgers’ equation with a+f3 sinx 
as “initial” condition, if y is identified as time t. Hence exact solution to the steady 
state can be obtained. This experiment simulates a 2D boundary layer problem 
because the shock “dissolves” near the boundary y = 0. We can adjust the thickness 
of the boundary layer (i.e., the smooth region near y = 0) by adjucting (3. We 
used ENO-RF-3 to compute (5.6), running it to steady state, with initial condition 
u(x,y,0) = 0 ! + /?sinx. We imposed periodic boundary conditions in x, enforced 
(5.6b) at y = 0, and imposed no boundary condition at the outflow boundary 
y = 1. (For a detailed discussion of th implementation of boundary conditions for 



- 20 - 


ENO schemes we refer the readers to [4]). Figures 29 and 30 contain the level curves 
for a = 0, 0 = 5 and a = 2, (3 = 5. 30 x 30 spatial points are used. □ 

EXAMPLE 5. This example was suggested by Professor R. Sanders (private 
communication), to study the smearing of contact discontinuities in different direc- 
tions when Ax / Ay. We solve 


Ut + U x + Uy = 0 , 



u is periodic in x and y with periods 1 
if (x-i) 2 + (y-i) 2 <! 
otherwise in 0 < x, y < 1 


(5.7a) 

(5.7b) 


with Ay = 2 Ax = and ran it for t = 2 (two periods). In the Figures 31-33 
the ‘+* signs are for the cross section y = \, and the signs are for the cross 
section x = We observe that ENO-Roe-3 does smear differently in the x and y 
directions with Ax ^ Ay, but ENO-Roe-A-3 does much better. On the other hand, 
ENO-Roe-S-2 (naive generalization to 2D) does not work well. It seems that some 
truly 2D subcell resolution techniques are needed. □ 


EXAMPLE 6. We consider the Riemann problems for the Euler equations of 
gas dynamics for a polytropic gas, i.e. (1.1) with d = 1, m = 3, and 


u = (p, M,E) t , f(u) = qn + (0, P,qP) T (5.8a) 


where 

P = -^pq 2 ), M = pq 


for the initial condition 



x < 0 
x > 0 


(5.8b) 


(5.8c) 


7 = 1.4 is used. For details of the Jacobian, its eigenvalues and eigenvectors, 
see [4], [11], or Example 9 below. 
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We use the same two sets of initial conditions as in [4]: 

{pl,Ql,Pl) = { 1,0,1); = 0.125,0,0.10) (5.9a) 

[pl^l.Pl) = (0.445,0.698,3.528); (5.9b) 

( Pr,Qr,Pr ) = (0.5,0,0.571) 


The results are in Figures 34-39. Notice the improved treatment of contact 
discontinuities after using the two contact sharpening algorithms of Section 3. Also 
notice that the comers of rarefaction waves (discontinuities in derivatives) are not 
resolved as well. Currently this phenomenon is still under investigation. □ 


EXAMPLE 7. This is the same equation as in Example 6 with the initial 


condition 


where 


u(z,0) 


Ul, 

' U M, 


0 < x < 0.1 
0.1 < x < 0.9 
0.9 < x < 1 


(5.10a) 


PL =PM = PR = 1 <1L = <]M = QR = 0, Pl = 10 3 , (5.10b) 

Pm = 10 - 2 , Pr = 10 2 . 

A solid wall boundary condition is applied to both ends. See [4], [14] for details. 


The results at the final time t = 0.038 are in Figures 40-45. The solid lines 
are the numerical solutions using ENO-RF-S-3 with 800 points. It can be regarded 
as an exact solution. Notice that ENO-RF-3 has essentially converged for the 
velocity and the pressure, but not for the density, due to the smearing of contact 
discontinuities. Also notice that ENO-RF-S-3 and ENO-RF-A-3 with 200 points 
have better resolution for the density than does ENO-RF-3 with 400 points. □ 
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EXAMPLE 8. In the above problems we only show pictures for third order 

schemes. If we compare with second and fourth order schemes we see the expected 

improvements in resolution. These improvements are usually not very significant 

for problems lacking structure in smooth regions, as is the case for most of the above 

examples. This is particularly true when the two contact-discontinuity-sharpening 

algorithms in Section 3 are used. However, the advantage of higher order schemes 

becomes significant, when the problem involved has some structure. To exemplify 

this, we solve the Euler’s equation in Example 6 with a moving Mach = 3 shock 

interacting with sine waves in density, i.e., initially 

f p = 3.857143; q = 2.629369; P = 10.33333 when x < -4 
< , (5.11) 

( p = 1 + esin5x; q = 0; P — 1 when x > — 4 

If e = 0, this is a pure Mach = 3 shock moving to the right. 

We take s = 0.1. For a linearized analysis of this problem see [8]. The results 
are in Figures 46-52. The solid lines are numerical solutions of ENO-RF-3 with 
1600 grid points. It can be regarded as the exact solution. We observe that the fine 
structure in the density profile makes the higher order schemes perform much better 
than the lower order methods. ENO-RF-3 with 200 points has roughly the same 
resolution as the second order MIJSCL type TVD scheme [9] with 800 points. On 
the other hand, the improvement of ENO-RF-3 over the second order TVD scheme 
is not so significant for the velocity and pressure profiles, because they both lack 
any detailed structure. □ 

EXAMPLE 9. We apply ENO schemes to 2D Euler’s equation of gas dynamics, 
i.e. (1.1) with d = 2, m = 4, and (we use f,g,x,y instead of f lf f 2 , xi,x 2 ): 

\l={p,M x ,M y ,E) T , f(u) = q x n + (0,P,0,q x P), (5.12a) 


g(u) = q y u + {0,0,P,q y P) 
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where 


p = (7 - i)(P - ^pg 2 ), g 2 = g 2 + gj, M x = pg s , M y = pg y 

Also 

c 2 = fr-i)(ir-i, 2 ). = ^ 

For §^, the eigenvalues are 


7i — gx c, A 2 — A3 — q x , A4 — gx "H c, 


the right eigenvectors are 


( 1 > 


(°\ 

g s-c 

, r 2 = 

0 

q y i 

1 

\H — q x c) 


My' 


and the left eigenvectors are 



( 1 ^ 


( 1 > 

q s ] 

, r 4 = 

q x + c 

Qy j 

1 Qy 

\\q 2 ) 


\H + q x c ) 


— biq x , 

C 


ti = (~gyiO,l,0), lz=[\-b2, b\q x , biq v , -61), 

£4 = r(^2 — biq x , —biq v ,bi) 


where: 


bi = 


7-1 


e “ 

62 = ^g 2 &i 


We can get the results for §& by symmetry. 


(5.12b) 

(5.12c) 

(5.13a) 

(5.13b) 

(5.13c) 

(5.13d) 

(5.13e) 


The test problem we choose is a moving shock interacting with compressible 
turbulence [17], [18]. At t = 0, a Mach 8 shock at x = —1.0 is moving to the right 
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into a state with Pr = 1, p R = 1 and q s = — ^-sin0.R cos(xkR cosd R + yk R sin##), 
q y = j^cosO R cos(xk R cosO R + yk R smO R ) where k R = 2tt, 9 r = §. We display 
the results at t = 0.20 in Figures 53-58. Notice that in [17], [18] similar results 
were obtained using a shock-fitting rather than a shock capturing method. This 
problem is difficult for shock capturing methods because both properties: high order 
accuracy (to resolve the fine structure) and an oscillation-free shock transition, are 
needed. A small-spurious oscillation will be enough to kill any resolution of the 
specified small fluctuations. This example shows that ENO schemes have excellent 
possibilities for shock-turbulence computations. □ 

VI. CONCLUDING REMARKS 

ENO schemes based on fluxes and TVD Runge-Kutta type time discretizations 
seem to work very well in our numerical tests which include ID and 2D scalar and 
systems problems. ENO-RF (Algorithm 2.3) seems preferred, since this method is 
half as expensive as ENO-LLF or ENO-LF, gives better resolution and seems always 
to converge to the correct entropy solution. To sharpen contact discontinuities, Al- 
gorithm 3.1 (for ID) or Algorithm 3.2 (for ID or 2D) can be used. High order 
ENO schemes show their special advantages for problems which have both discon- 
tinuities and detailed fine structure in smooth regions e.g. Example 8 and Example 
9 in Section V. Among the practical things needing further investigation are the 
resolution of rarefaction corners (discontinuities in derivatives). Finally theoretical 
justification for the evident stability of these methods would be quite welcome. 
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Conceming the figures: 

(1) In Figures 1-18 and 31-39, the solid lines are the exact solutions, and the circles 
or diamonds are the numerical solutions. 

(2) Figures 19-28 are the surfaces of the numerical solutions. 

(3) Figures 29-39 and 53-58 are the level curves of the numerical solutions. 

(4) In Figures 40-45, the solid lines are the numerical solution of ENO-RF-S-3 
with 800 points, and the diamonds are the numerical solutions explained in the 
legends. 

(5) In Figures 46-52, the solid lines are the numerical solution of ENO-RF-3 with 
1600 points, and the diamonds are the numerical solutions explained in the 
legends. 



Figure 1: 
Figure 2: 
Figure 3: 

Figure 4: 
Figure 5: 
Figure 6: 
Figure 7: 
Figure 8: 
Figure 9: 
Figure 10: 
Figure 11: 
Figure 12: 
Figure 13: 
Figure 14: 
Figure 15: 
Figure 16: 
Figure 17: 
Figure 18: 
Figure 19: 

Figure 20: 

Figure 21: 
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Figure legends 

ENO-RF-3, (5.1) with ui e ft = 2, u r ; g ht = -2 
ENO-RF-3, (5.1) with t*i e ft = -3, unght = 3 
ENO-RF-3, (5.2) with u°(x) = 1 in [ — 0], and 

u°(x) = 0 elsewhere 
ENO-Roe-3, (5.3)-(5.4a) 

ENO-Roe-S-3, (5.3)-(5.4a) 

ENO-Roe-A-3, (5.3)-(5.4a) 

ENO-Roe-3, (5.3)-(5.4b) 

ENO-Roe-S-3, (5. 3) -(5. 4b) 

ENO-Roe-A-3, (5.3)-(5.4b) 

ENO-Roe-3, (5.3)-(5.4c) 

ENO-Roe-S-3, (5.3)-(5.4c) 

ENO-Roe-A-3, (5.3)-(5.4c) 

ENO-Roe-3, (5.3)-(5.4d), t = 2 
ENO-Roe-S-3, (5.3)-(5.4d), t = 2 
ENO-Roe-A-3, (5.3)-(5.4d), t = 2 
ENO-Roe-3, (5.3)-(5.4d), t = 8 
ENO-Roe-S-3, (5.3)-(5.4d), t = 8 
ENO-Roe-A-3, (5.3)-(5.4d), t = 8 

ENO-RF-3, (5.5) with ( 111 , 1 * 2 , 1 * 3 ,^ 4 ) = (—1,— 0 . 2 , 0.5, 0 . 8 ), 
t = 1 , 20 x 20 points 

ENO-RF-3, (5.5) with (ui,u 2 ,u 3 ,u 4 ) = (- 1 , - 0 . 2 , 0.5, 0 . 8 ), 
t = 1 , 80 x 80 points 

First order Engquist-Osher scheme, (5.5) with 

(t*i, t* 2 , U 3 , u 4 ) = ( — 1, — 0.2, 0.5, 0.8) , t — 1, 80 x 80 points 
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Figure 22: 

Figure 23: 

Figure 24: 

Figure 25: 

Figure 26: 

Figure 27: 

Figure 28: 

Figure 29: 
Figure 30: 
Figure 31: 
Figure 32: 
Figure 33: 
Figure 34: 
Figure 35: 
Figure 36: 
Figure 37: 
Figure 38: 
Figure 39: 
Figure 40: 
Figure 41: 


ENO-RF-3, (5.5) with (ux,tt 2 ,U 3 ,u 4 ) = (-0.2, -1,0.5, 0.8), 
t = 1, 80 x 80 points 

ENO-RF-3, (5.5) with (tii,u 2 ,« 3 ,U 4 ) = (0.5,— 1,— 0.2, 0.8), 
t = 1, 80 x 80 points 

ENO-RF-3, (5.5) with (ui,ti 2 ,U 3 ,U 4 ) = (—1,0.5, —0.2, 0.8), 
t — 1 , 80 x 80 points 

ENO-RF-3, (5.5) with (ux,ti 2 >“ 3 >« 4 ) = (—1,— 0.2, 0.8, 0.5), 
t = 1, 80 x 80 points 

ENO-RF-3, (5.5) with (tix,U 2 >« 3 >« 4 ) = (0.8, —1, —0.2, 0.5), 
t = 1, 80 x 80 points 

ENO-RF-3, (5.5) with (ux,U 2 ,U 3 ,« 4 ) = (0.8, —1,0.5, —0.2), 
t = 1, 80 x 80 points 

ENO-RF-3, (5.5) with («x,« 2 > « 3 >“ 4 ) = (0.8, —0.2, —1, 0.5), 
t = 1 80 x 80 points 

ENO-RF-3, (5.6) with a = 0, (3 = 5, 30 x 30 points 

ENO-RF-3, (5.6) with a = 2, 0 = 5, 30 x 30 points 

ENO-Roe-3, (5.7), '+’ for y = i, for x = £, t = 2 
ENO-Roe-A-3, (5.7), '+’ for y = |, for x = |, t = 2 
ENO-Roe-S-2, (5.7), V for y = i, '** for x = t = 2 
ENO-RF-3, (5.8)-(5.9a), 100 grid points, t — 2, density 
ENO-RF-3, (5.8)-(5.9a), 100 grid points, t = 2, velocity 
ENO-RF-3, (5.8)-(5.9a), 100 grid points, t = 2, pressure 
ENO-RF-3, (5.8)-(5.9b), 100 grid points, t = 1.3, density 
ENO-RF-S-3, (5.8)-(5.9b), 100 grid points, t = 1.3, density 
ENO-RF-A-3, (5.8)-(5.9b), 100 grid points, t = 1.3, density 
ENO-RF-3, (5.8a,b)-(5.10), 200 points t = 0.038, density 
ENO-RF-3, (5.8a,b)-(5.10), 400 points, t = 0.038, density 
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Figure 42: 
Figure 43: 
Figure 44: 
Figure 45: 
Figure 46: 
Figure 47: 
Figure 48: 

Figure 49: 
Figure 50: 

Figure 51: 
Figure 52: 

Figure 53: 
Figure 54: 
Figure 55: 
Figure 56: 
Figure 57: 
Figure 58: 


ENO-RF-3, (5.8a,b)-(5.10), 400 points, t = 0.038, velocity 
ENO-RF-3, (5.8a,b)-(5.10), 400 points, t = 0.038, pressure 
ENO-RF-S-3, (5.8a,b)-(5.10), 200 points, t = 0.038, density 
ENO-RF-A-3, (5.8a,b)-(5.10), 200 points, t = 0.038, density 
ENO-RF-3, (5.8a,b)-(5.11), 200 points, t = 1.8, density 
ENO-RF-3, (5.8a,b)-(5.11), 400 points, t = 1.8, density 
Second Order MUSCL Type TVD Scheme, (5.8a,b)-(5.11), 
800 points, t = 1.8, density 

ENO-RF-3, (5.8a,b)-(5.11), 200 points, t = 1.8, velocity 
Second Order MUSCL Type TVD Scheme, (5.8a,b)-(5.11), 
200 points, t = 1.8, velocity 

ENO-RF-3, (5.8a,b)-(5.11), 200 points, t = 1.8, pressure 
Second Order MUSCL Type TVD Scheme, (5.8a,b)-(5.11), 
200 points t = 1.8, pressure 

ENO-RF-3, Example 9, 60 x 40 points, t = 0.2, pressure 
ENO-RF-3, Example 9, 60 x 40 points, t = 0.2, vorticity 
ENO-RF-3, Example 9, 60 x 40 points, t = 0.2, entropy 
ENO-RF-3, Example 9, 120 x 80 points, t = 0.2, pressure 
ENO-RF-3, Example 9, 120 x 80 points, t = 0.2, vorticity 
ENO-RF-3, Example 9, 120 x 80 points, t = 0.2, entropy 
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Figure 20: ENO-RF-3 , (5.5) with (u. ,u-,u-,u / )»(-l,-0.2,0.5,0.8) , 

t = 1, 80 x 80 points 
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Figure 21: 
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First Order Engqulst-Osher Scheme, (5.5) with 
(u^ f U2>u^ > u^)~(-l f -0.2 t 0.5 > 0.8), t-1, 80 x 80 points 
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Figure 37: ENO-RF-3, (5.8)— (5. 9b) , 100 points, t"1.3, density 



Figure 38: ENO-RF-S-3, (5.8)-(5.9b) , 100 points, t”1.3, density 
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Figure 42: EN0-RF-3, (5.8a,b)-(5. 10) , 400 points, t-0.038, velocity 
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Figure 46: EN0-RF-3, (5.8a,b)-(5. 11) , 200 points, . t-1. 8, density 


4 




-4 -2 ' 0 2 4 

Figure 48: Second Order MUSCL Type TVD Scheme, (5.8a,b)-(5. 11) , 800 points 

t=1.8, density 
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-4 - 2 ' 0 2 
Figure 49: ENO-RF-3, * (5.8a,b)-(5. 11) , 200 points, t-1.8, velocity 
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Figure 52: Second Order MUSCL Type TVD Sch.jie, (5. 8a, b)-(5. 11) , 200 points, 

t=1.8, pressure 







0.97 A 1 ' I • ‘ ■ i > ' ■ ' l 1 

*1.47 -0.74 0.00 0.74 1.47 

Figure 53: ENO-RF-3, Example 9, 60 x 40 points, t«0.2, pressure 



Figure 54: ENO-RF-3, Example 9, 60 x 40 points, t«0.2, vorticity 
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Figure 55: ENO-RF-3, Example 9, 60 x 40 points, t m 0.2, entropy 
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Figure 56: EN0-RF-3, Example 9, 120 x 80 points, t»0.2, pressure 
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Figure 58: ENO-RF-3, Example 9, 120 x 80 points, 't=0.2, entropy 
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